function [g gt gtt] = grav_metric(params, u, udt, udtdt)



  % load precomputed function for order = 2 computation
  if params.metricOrder == 2
    if params.metricLinearFnInterpolation
      global nh2s7 nk2s7 na2t2 na2t4 na2t5;
      global interpgrid;

      load('pde/nflin.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5', 'z');
      interpgrid = z;

    else
      global nh2s7 nk2s7 na2t2 na2t4 na2t5;
      if params.metricMark
        load('pde/nfmark.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');
      else
        load('pde/nf.mat',     'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');
      end
    end
  end
  % assumptions: all of the precomputed functions are on a chebyshev grid, from (0, 2). 


  Nt = size(u, 1);
  tord = Nt - 1;
  tord


  % **********************
  %  compute g
  % **********************

  disp(sprintf('computing metric at order %d', params.metricOrder));

  % allocate memory for intermediate g computation:
  g = zeros(Nt, params.Nx, params.Ny, params.Nz, 4, 4);
  

  for ti=1:Nt
    tic;
    g(ti, :, :, :, :, :) = grav_metric_eval(u(ti,:), udt(ti,:), udtdt(ti,:), params);
    gtime = toc;
    disp(sprintf('computed metric at t=%f in %.1fsec', params.t, gtime));
  end


  % **********************
  %  compute time derivatives of g
  % **********************

  gt  = zeros(params.Nx, params.Ny, params.Nz, 4, 4);
  gtt = zeros(params.Nx, params.Ny, params.Nz, 4, 4);

  d1k = d1ker(tord);
  d2k = d2ker(tord);

  gt  = conv2(g(:, :), d1k, 'valid')./params.dt;
  gtt = conv2(g(:, :), d2k, 'valid')./params.dt./params.dt;

  % ******************************
  %  splice out intermediate g
  % ******************************

  g = g(tord/2+1, :, :, :, :, :);
  gt = reshape(gt, 1, params.Nx, params.Ny, params.Nz, 4, 4);
  gtt = reshape(gtt, 1, params.Nx, params.Ny, params.Nz, 4, 4);
end



function ker = d1ker(order)
  switch (order)
    case 2
      ker = [-1/2; 0; 1/2];
    case 4
      ker = [1/12; -2/3; 0; 2/3; -1/12];
    case 6
      ker = [-1/60; 3/20; -3/4; 0; 3/4; -3/20; 1/60];
    case 8
      ker = [1/280; -4/105; 1/5; -4/5; 0; 4/5; -1/5; 4/105; -1/280];
  end
  ker = flipud(ker);
end


function ker = d2ker(order)
  switch (order)
    case 2
      ker = [1; -2; 1];
    case 4
      ker = [-1/12; 4/3; -5/2; 4/3; -1/12];
    case 6
      ker = [1/90; -3/20; 3/2; -49/18; 3/2; -3/20; 1/90];
    case 8
      ker = [-1/560; 8/315; -1/5; 8/5; -205/72; 8/5; -1/5; 8/315; -1/560];
  end
end

function udtdt = exactdtdt(params, hydro_sol, ts)

  udtdt = zeros(params.Nx*params.Ny*5, length(ts));

  timestep = params.Lx * 0.0000001;

  tord = params.timeOrder;  % order of time derivatives 
  d1k = d1ker(tord);
  for ti=1:length(ts)
    ts2 = [-tord/2:1:tord/2] .* timestep + ts(ti);
    [u, udt] = deval(hydro_sol, ts2);
    tmp = conv2(udt, d1k', 'valid');
    udtdt(:, ti) = conv2(udt, d1k', 'valid');
  end
  udtdt = udtdt./timestep;
end

